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A large fraction of cosmological information on dark energy and gravity is encoded in the nonlinear 
regime. Precision cosmology thus requires precision modeling of nonlinearities in general dark energy 
and modified gravity models. We modify the Gadget-2 code and run a series of N-body simulations 
on modified gravity cosmology to study the nonlinearities. The modified gravity model that we 
investigate in the present paper is characterized by a single parameter ^, which determines the 
enhancement of particle acceleration with respect to general relativity (GR), given the identical 
mass distribution ((^ = 1 in GR). The first nonlinear statistics we investigate is the nonlinear 
matter power spectrum at fc < 3/i/Mpc, which is the relevant range for robust weak lensing power 
spectrum modeling at i" < 2000. In this study, we focus on the relative difference in the nonlinear 
power spectra at corresponding redshifts where different gravity models have the same linear power 
spectra. This particular statistics highlights the imprint of modified gravity in the nonlinear regime 
and the importance of including the nonlinear regime in testing GR. By design, it is less susceptible 
to the sample variance and numerical artifacts. We adopt a mass assignment method based on 
wavelet to improve the power spectrum measurement. We run a series of tests to determine the 
suitable simulation specifications (particle number, box size and initial redshift). We find that, the 
nonlinear power spectra can differ by ~ 30% for 10% deviation from GR = 0.1) where the rms 

density fiuctuations reach 10. This large difference, on one hand, shows the richness of information 
on gravity in the corresponding scales, and on the other hand, invalidates simple extrapolations of 
some existing fitting formulae to modified gravity cosmology. 

PACS numbers: 98.65.Dx,95.36.+x,04.50.+h 



I. INTRODUCTION 

One of the biggest challenges of modern cosmology and 
physics is the existence of the dark universe. Assuming 
the validity of general relativity (GR), cosmological ob- 
servations lead to the discovery of dark matter and dark 
energy, which account for ^ 96% of the total matter and 
energy budget of the Universe (e.g. [l[). However, since 
we do not have independent tests of GR at relevant scales, 
the same set of observations could imply another pos- 
sibility, the failure of general relativity at galactic and 
cosmological scales. This possibility, which serves as an 
alternative to dark matter/dark energy, has become an 
area of active research. Discriminating between the dark 
matter/dark energy and modified gravity (MG) models, 
testing GR at cosmological scales and probing dark mat- 
ter and dark energy through cosmological observations, 
are thus an entangled task, of crucial importance for both 
cosmology and physics. 

Challenges exist in both the observation side and the- 
ory side. Although there are numerous and potentially 
powerful observations suitable for this task Q , their 
precision measurements are challenging. On the other 
hand, much of the cosmological information is encoded 
in the nonlinear regime. Modeling the nonlinearities to 
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the required ^ 1% accuracy is challenging too, even for 
the simplest case, the standard ACDM cosmology with 
only gravitational interaction (e.g. 0-0|)- Cosmologies 
based on dynamical dark energy or MG are facing similar 
requirements. References @t1I1 have performed N-body 
simulations for dynamical and coupled dark energy mod- 
els [12]. 

Comparing to the dark matter/dark energy cosmology, 
understanding the evolution of the Universe in MG mod- 
els is often more difficult, due to the intrinsically nonlin- 
ear feature of gravity in these models or the existence of 
extra dynamical fields. Despite these difficulties, the ex- 
pansion history of the Universe and the structure growth 
to the first order have been robustly understood for many 
of the MG models such as TeVeS JliJ3, DGP(short 
for Dvali,Gabadadze and Porrati) |l5l, 16| and the /(i?) 
gravity |17H2lj . People have also achieved success in 
understanding the nonlinear evolution through analyti- 
cal and semianalytical methods (e.g. 122-25]) Recently, 
self-consistent gravity solvers for f {R) [26l - l29| and DGP 
gravity [3Ci| models have been developed and led to sig- 
nificantly improved understanding of the nonlinear evolu- 
tion. Simulations with extra scalar fields and interaction 
with dark matter have also been performed (e.g. [sij). 

Since deviations from GR in general lead to nonlinear 
differential equations of gravity, in principle we have to 
develop the suitable N-body codes for each viable MG 
model and run the corresponding simulations. However, 
since we do not have the final theory of gravity based 
on the first principles, there are in principle infinite MG 
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models to be investigated. One possibility to circumvent 
this problem is to choose a suitable parameterization for 
the MG models and run a finite number of simulations 
to sample the relevant parameter space. We are then 
able to interpolate/extrapolate the simulation results to 
explore the whole relevant parameter space. 

The statistics we focus on is the matter power spec- 
trum. It determines the lensing power spectrum and is 
also highly relevant to the 2D galaxy clustering, both 
will be measured to high precision by ongoing and planed 
imaging surveys, such as DES, LSST, JDEM/SNAP, Eu- 
chd/DUNE, and KDUST. As we will show later, the 
density evolution is determined by a single parameter 
of gravity, C, which quantifies the ability of mass con- 
centration to distort the space-time metric. In princi- 
ple, C can be both scale, time, and environmental de- 
pendence. Comprehensive investigation on general C is 
beyond the scope of the current paper. Instead, we will 
adopt a highly simplified form of C and run a series of 
N-body simulations to quantify the nonlinear evolution 
of the Universe. 

As shown in Heitmann et al. 2008, (hereafter H08, 
Q), however, to run simulations and model nonlinear 
matter power spectrum to 1% accuracy to k ^ 1/iMpc 
is very challenging, requiring Gpc or larger simulation 
box size, 1024^ or more particles, and beyond. Being 
aware of these difficulties and limited computation re- 
source, we take a modest goal, to quantify the influence 
of MG on the nonlinear matter power spectrum with re- 
spect to the standard ACDM to ~ 1% accuracy. Namely, 
the statistics that we will focus on is the relative differ- 
ence between the nonlinear matter power spectra in the 
given MG model and in ACDM. We will choose the right 
redshifts of simulation output such that the linear matter 
power spectra in the given MG models are equal to the 
ones in the corresponding ACDM. This particular statis- 
tics has a number of attractive features. First, it isolates 
and highlights the role of MG in nonlinear evolution. Sec- 
ond, it reduces much of the numerical artifacts by taking 
ratios. Similar tricks have been adopted in many pre- 
vious simulations (e.g. Third, it can improve 
the efficiency to understand nonlinearities in MG mod- 
els, which can now be reduced to two separate ingredi- 
ents: the nonlinear evolution in ACDM and the relative 
difference between MG models ACDM. 

The current paper only analyzes a very limited set of 
simulations. Nonetheless, it robustly show that, even 
after scaling out the difference in the linear evolution, 
gravity still leaves significant features in the nonlinear 
power spectrum. In subsequent studies, we will run sim- 
ulations covering larger parameter space to better un- 
derstand these features and hopefully develop a general 
fitting formula. Furthermore, we will study the peculiar 
velocity power spectrum and the redshift distortion (3D 
galaxy clustering), based on these simulations. The ongo- 
ing spectroscopic redshift surveys like BOSS, LAMOST 
and WiggleZ, and planned spectroscopic redshift surveys 
like BigBOSS, JDEM/ADEPT, Euclid/SPACE and SKA 



are able to measure these statistics to unprecedented ac- 
curacy. We will also investigate the halo statistics, one 
of the key scientific goals for galaxy and cluster surveys. 

This paper is organized as follows. In Sec. |lll we 
present the MG parameterization adopted for the simu- 
lations, the precision requirements and the code specifi- 
cations. In Sec. mil we test the accuracy of the simula- 
tions. We present major simulation results in Sec. IIVI , 
discussion in Sec. IVland more results in the Appendix. 



II. THE SIMULATION LAYOUT 
A. The (J parameterization on modified gravity 

There are several existing parameterizations and gen- 
eral guidances on modified gravity [1, HI, Isl - fssj . What 
we adopt in this paper is the C, parameterization. It is a 
condensed version of the GcS-rj parametrization 0, [s^ , 
which quantifies two key aspects of gravity. 

To understand this point, we begin with the structure 
formation in GR, for which the central issue is to deter- 
mine the particle acceleration given the mass distribu- 
tion. The scalar perturbation of the space-time metric is 
described by two potentials, ds^ = — (1 + 2'ip)dt^ -\- a^(l -I- 
20)c?x^. The usual Poisson equation, k'^(t) — AirGa^Sp 
(in Fourier space), relates the potential (j) to the matter 
distribution p. However, (j) is not the potential directly 
responsible for the structure formation in our N-body 
simulations of nonrelativistic cold dark matter particles. 
The contribution to the particle acceleration from this 
potential is suppressed by a factor (f/c)^ ^ 1, compar- 
ing to the contribution from the other potential iJj. Thus 
for the nonrelativistic cold dark matter particles that our 
simulations deal with, their acceleration is determined 
solely by -07 d{av)/dt = ikV', where v is the proper mo- 
tion. GR predicts ip — —cf), if dark energy anisotropic 
stress is negligible. Now, given an initial mass distribu- 
tion p, we obtain cj) from the Poisson equation, with the 
coupling constant G. Then through the relation tp = 
we obtain ip and then the acceleration. Thus given the 
initial positions (density) and velocities of particles, we 
can move particles in each simulation time step and then 
have a closed procedure to simulate the evolution of the 
Universe under gravity. 

A natural parametrization of modified gravity is thus 
to replace the Newton's constant G by the effective 
Newton's constant GeS and the relation ip — —0 by 
77 = —(j}/ip. Now, given the mass distribution, the ac- 
celeration is solely determined by the combination ( [t^I , 
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r]{k,z) 

This is the quantity that enters into the tp-p relation, 

k^iP = -CAnGdp . (2) 

GR has the value C = 1- Clearly, if the two universes 
have the identical initial conditions, identical expansion 
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rate and identical C(fc, z), the statistics of the density and 
velocity fields would be identical [73]. 

Thus, instead of running a series of simulations on 
a 2D grid of GcS-r] parameter space, we just need to 
run a series of simulations on a ID grid of ^ parameter 
space. It significantly reduces the amount of simulations 
required. This is the major reason that we adopt this sin- 
gle parameter parametrization on modified gravity. Be- 
sides it, there are a number of attractive features of this 
parametrization. 

First of all, many MG models, such as the DGPgrav- 
ity model il_6i] and the f{R) gravity model [3 [31 (i^ 
the linear regime), the Yukawa-like MG model and the 
7-index MG model j33| fit into this parameterization. 
Second, such parameterization requires minimum mod- 
ification in the N-body gravity solver, does not require 
extra computation time, and thus is suitable for fast ex- 
ploration of the vast parameter space of MG models. In 
fact, there already exists a number of simulations on 
Yukawa- like gravity [s^ Third, Gcff and 77 (and 

hence C) , can be measured in a rather model independent 
manner, by combining imaging surveys and spectroscopic 
surveys [la, EH . This links theories and observations di- 
rectly. Furthermore, the reconstruction accuracy can be 
improved by including all available data and performing 
a multiparameter fitting [H, . 

Clearly, this parameterization does not capture all fea- 
tures of MG, such as the environmental dependence of 
gravity, as found in f{R) gravity [44| and the DGP model 
|45| . Nevertheless, the simulations based on this param- 
eterization serve as an useful step toward better under- 
standing of MG cosmology. The simulation results can 
be used as templates to understand more complicated 
MG models. A close analogy is the scale-free simula- 
tions. Although the real CDM(cold dark matter) trans- 
fer function is certainly not scale- free (power-law), these 
scale- free simulations do significantly improve our under- 
standing of the nonlinear evolution of structure forma- 
tion. They are helpful in developing fitting formula like 
that of Peacock-Dodds (hereafter PD96, [4g|) and Smith 
et al. 2003 (hereafter halofit, [ij). We hope that sim- 
ilar procedure applies to the case of MG models. For 
example, the formalism proposed by [l^ relies on the 
interpolation between the nonlinear power spectrum in 
GR and the one in MG without environmental depen- 
dence. Understanding the nonlinearities in MG models 
without environmental dependence thus serves as a natu- 
ral step to understand nonlinearities in more complicated 
MG models. 

Modifying existing N-body codes to incorporate the ( 
parameterization is straightforward. The only modifica- 
tion is to change the particle acceleration a to C x a- In the 
simulation setup, we fix the expansion rate identical to 
that of the flat ACDM cosmology [tI] . In addition, we do 
not aim to explore the whole space of C(fc, z). Rather, we 
will focus on very special cases of C and postpone the gen- 
eral investigation for future studies. The z) adopted 
in our simulations is scale independent (C(fc,z) — C{z))- 



The success of CMB (cosmic microwave background) and 
BBN (big-bang nucleosynthesis) implies that GR is likely 
valid in the early Universe. For this reason, we adopt a 
step function in z, such that ^ = 1 at z > zmg and 
C =constant7^ 1 at z < zmg- Throughout this paper, we 
have adopted zmg = Zi — 100, where z^ is the initial red- 
shift of simulations. Since we have GR valid at high red- 
shift (z > 100), the transfer function at Zi = 100 adopted 
in the MG models is identical to that in GR. For the 
adopted MG parameterization, the linear density growth 
factor D{k,z) is scale independent D{k,z) = D{z). Thus 
the linear power spectrum for modified gravity models 
only differs from ACDM by the linear density growth 
factor D{z). In a companion paper, we will explore MG 
models with other redshift dependence. 



B. The precision requirements 

All MG simulations begin with the identical initial con- 
dition at Zi = 100. Since the adopted C is scale indepen- 
dent, the linear density growth factor D{z^Q is scale in- 
dependent, as can be seen from the equation at z < zmgj 
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Here, = d/da and = d^/da?'. fio is the present day 
matter density in unit of the critical density, i/o and 
H are the present day Hubble constant and the Hubble 
parameter at z = 1/a — 1. (5m is the linearly evolved 
matter over-density and D cx 5m is the linear density 
growth factor. Thus, given a redshift zg in the standard 
ACDM, we can find the corresponding redshift Z(^ in the 
MG universe, such that 



D{zsX^l) = Diz^X) 



(4) 



Here, the subscript S denotes the standard ACDM cos- 
mology. Since all the simulations begin with the identical 
initial condition, the above relation means that, 

PL(fc;z5,C = 1) = PL{k;zc,0 ■ 

Here Pl is the linear matter power spectrum. Through- 
out this paper, we use the subscript "L" for the linear 
statistics and the subscript "NL" for the nonlinear statis- 
tics. 

Modifications in GR change the structure growth his- 
tory. The structure grows faster in a universe with big- 
ger C.. The primary quantity that we want to measure 
through the simulations is 



^ PNL(fc;^:c^C) 

" PNL(fc;2:s,C = l) 



(5) 



e ^ 1 has a number of implications. (1) If the nonlinear 
power spectrum is completely determined by the linear 
one, independent of the expansion and structure growth 



4 



history and the underlying gravity, then e = 1 . A number 
of fitting formulae applicable to GR have been extended 
to study the nonlinear evolution in MG models, based 
on this assumption. Thus e provides a direct test on the 
applicability of these fitting formulae to MG models. Pre- 
cision cosmology requires that, only if|e — l|<10^^in 
the relevant k range, may the systematical error induced 
by these fitting formulae be subdominant. Otherwise, 
significant modifications shall be made. (2) e / 1 also 
means that there is extra information of gravity encoded 
in the nonlinear matter power spectrum, which does not 
show up in the linear power spectrum at the same epoch. 
This helps to test GR at nonlinear regimes. Such in- 
formation is complementary to those in the linear power 
spectrum at the same epoch and those in the deeply non- 
linear regime where gravity reduces to GR through en- 
vironmental dependence mechanisms like the chameleon 
mechanism and the Wainshtein mechanism [isj . 

Much of the cosmological information in weak lens- 
ing surveys come from the lensing power spectrum mea- 
surement at < 2000 of source galaxies at ~ 1. 
Since the lensing kernel peaks at half way between the 
source and the observer, the peak contribution comes 
from k ~ e/[x{zs)/2] < 2/i/Mpc. At £ = 2000, the sta- 
tistical error in the lensing power spectrum measurement 
can reach below 1% for the planning of wide surveys. Un- 
der the Limber approximation, the lensing angular (2D) 
power spectrum is linearly proportional to the 3D nonlin- 
ear matter power spectrum. Thus, to match the observa- 
tion accuracy, we set a goal to model e to ~ 1% accuracy 
at z ~ 0.5 and k < 3/i/Mpc. 

Since the simulations run from the identical initial con- 
dition, the cosmic variances in the resulting power spec- 
tra Pjvi of different MG models are highly (positively) 
correlated. Since the simulations are run by the same 
code, with the same time steps, errors induced by the 
numerical artifacts into Pnl should also be highly (pos- 
itively) correlated. When taking the ratio of two power 
spectra to evaluate e, much of the errors in Pnl cancels. 
We thus expect higher accuracy in e than in Pnl- Thus, 
once we control the error in P^l to ^ 1% accuracy, we 
are likely able to measure e to 1% accuracy. 

We run a set of = 512"^ particle N-body simula- 
tions using the GADGET-2 code, on the 32-CPU Ita- 
nium server at the Shanghai astronomical observatory. 
All the simulations that we use to calculate e adopt 
L = 300h^^ Mpc. Adopting a smaller box size allows 
us to go deeper into the nonlinear regime. However, a 
smaller box size can cause numerical artifacts, due to the 
missing of power at k < 2tt/L, which affects the nonlin- 
ear evolution through mode coupling (5| . Another reason 
that we do not adopt a smaller box size is that, we plan to 
use the same simulations for velocity and halo statistics, 
which prefer a larger box size. 



C. The GADGET-2 simulation specifications 

We adopt a parallel GADGET-2 N-body code [H ^ 
to run the simulations. With a TreePM algorithm, where 
only short-range forces are computed with the "tree" 
method while long-range forces are determined by par- 
ticle mesh (PM) algorithm, GADGET-2 combines high 
efficiency with high resolution. 

The background expansion history is fixed as the one 
in a flat ACDM cosmology with the matter density 
Hq — 0.276 and the cosmological constant ft a = 0.724. 
The transfer function is fixed by the above parameters, 
the baryon density fif, = 0.046 and the dimensionless 
Hubble constant h = 0.703. The amplitude of the ini- 
tial fiuctuations is fixed such that, if linearly evolved to 
z = in the adoption of ACDM cosmology, the rms den- 
sity fluctuation within a sphere of radius 8/i~^Mpc is 
as = 0.811. 

We use 512'^ PM mesh grids through all the simula- 
tions. The force softening length 7 depends on the mean 
inter particle separation, with 7 = 0.022L/A^'^/'^, where 
L is the box size and N is the particle number. For simu- 
lations performed with 512'^ particles in the 300 h~^Mpc 
box, 7 = 12.89 /i-^kpc. In GAD GET-2, the adaptive 
time step is set by At = ^2^7/|a|, where ^ controls time 
step accuracy and a is the acceleration. ^ is fixed at 0.5% 
for all the simulations. With the adopted small softening 
length, the number of total adaptive time steps for our 
ACDM simulation is about 4000. Fig. 13 of H08 shows 
that, for 3000 time steps in total, the resulting difference 
in the power spectra is less than 0.04%. We thus believe 
that, the time stepping we adopt suffices for the purpose 
of this paper. 

III. SIMULATION TESTS 

In this section, we present steps to control the robust- 
ness of simulation results. We adopt the Daubechies mass 
assignment method to improve the accuracy of power 
spectrum measurement. We run a number of tests to 
justify that the adopted simulation specifications (parti- 
cle number, simulation box size and initial redshift) are 
adequate to constrain the nonlinear power spectrum out 
to fc = 3/iMpc~^ with - 1% accuracy. Finally, we show 
that the modified GADGET-2 code reproduces the cor- 
rect linear evolution in the linear regime. 

A. Calculating the matter power spectrum 

Usually people use the fast Fourier transform (FFT) 
to calculate the matter power spectrum. This requires 
assigning simulation particles to uniform grids first. For 
commonly used mass assignment methods, the resulting 
power spectrum is biased by the smoothing and aliasing 
effects, even at scales well below the Nyquist frequency 
(e.g. [sij). To reach the required accuracy, we must 
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FIG. 1: The top panel shows the matter power spectra for 
simulations of different particle numbers. The bottom panel 
shows the ratios with respect to 512^ particle simulation. 
All the power spectra are calculated through FFT on 512"' 
meshes. The dotted vertical line shows the scale of 0.7fcNy 
for our FFT power spectrum measurement. The black long 
dashed line is the linear power spectrum. At z = 0, nonlinear 
correction becomes significant at fc > 0.2/iMpc~^. 



FIG. 2; The impact of initial redshift on the matter power 
spectrum. In the top panel, we show the power spectra of 
the two simulations with starting redshift Zi — 49( red line) 
and Zi = 100 (the black line). The bottom panel shows the 
relative difference. The black horizontal dotted line shows 
the 1% precision requirement. The red vertical dotted line is 
0.7fcNy. The same as Fig. [1] the black long dashed line shows 
the linear power spectrum. 



correct for these biases. Reference [5l| proposes an iter- 
ative method to perform such task. Alternatively, [s^ 
adopts the Daubechies vifavelet transformation for the 
mass assignment. The scale function of the Daubechies 
wavelets transform has compact top-hat like support in 
the Fourier space, which avoids the sampling effect and 
allows computationally efficient mass assignment onto 
grids. Using this scale function to do the mass assignment 
allows for robust measurement of the power spectrum to 
52| . Throughout this paper, we will adopt 



k = 0.7k 



Ny 



this method to calculate the matter power spectrum. 



B. Particle number 

The particle number in GADGET-2 controls the mass 
resolution, force resolution and the time step. A larger 
particle number is necessary to avoid errors from dis- 
creteness effects at small scales of interest. As pointed 
out by Sirko fsS^, although simulations can probe the 
evolution of structures beyond the particle Nyquist fre- 
quency, fcNy.p — ttN^^^/L, it is unclear whether or not 
the shot noise term beyond this frequency already in the 
initial condition will impact power at the wavenumbers 
of interest. The issue may be made moot merely by us- 
ing negligible values of V/N in simulations. How many 
particles are required to sufficiently sample the density 
field and calculate the matter power spectrum robustly 
to A: = 3/i/Mpc? To answer this question, we run three 
simulations with identical initial conditions and a box 
size of aOO/i-^Mpc, but with 128^, 256^ and 512^ par- 
ticles, respectively. Fig. [1] shows the nonlinear power 



spectra calculated by the Daubechies' mass assignment 
method. We see clearly the impact of particle number on 
the simulated power spectrum in the nonlinear regime. 
The relative difference between the 256^ and 128^ results 
at l/i/Mpc< k < 3ft./Mpc is 4%, implying a minimum 
error of 4% in the 128"^ particle simulation, due to the 
resolution limitation. But the relative difference reduces 
to below 1-2% between the 512^ and 256^ ones, showing 
that the resolution induced error in the 256"^ particle sim- 
ulation is reduced significantly. This trend of convergence 
implies that the resolution induced error in the 512'^ sim- 
ulation is likely below ~ 1%. We then speculate that, if 
the Daubechies' mass assignment method was adopted, 
nonlinear power spectrum in the 512'^ particle simulation 
can attain 0(1%) accuracy out to fc ^ 3/iMpc~^. To ro- 
bustly test it, higher resolution simulations (e.g. ones 
with 1024'^ particles or more) are required. This test 
shall be performed in future works. 



C. Initial redshift 

Testing the effect of changing the starting redshift in 
simulations is also important. Since the initial condition 
is generated under the Zel'dovich approximation [s^l , the 
initial redshift can not be too low, otherwise higher or- 
der corrections can be non-negligible. However, it is not 
automatically the case that higher Zi is better, because 
numerical errors (most obviously suppression of power by 
limited force resolution) have more time to accumulate 
in that case Q. Our initial redshift tests are started at 
Zi = 49 and Zi = 100 respectively, both with 512'^ par- 
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FIG. 3: The linear density growth, normalized a.t z = Zi — 
100. The squares are obtained from the power spectra at 
k = 0.025 hMpc~^ in the simulations. The corresponding 
redshifts are shown in Table H] The difTerent colored and style 
lines are the theoretical predictions. The agreement between 
the simulated results and theoretical predictions justifies our 
modified GADGET-2 code for MG models. 

tides and a 300/i~^Mpc box size. Fig. [2] compares the 
two power spectra at redshift z — 0. The agreement is 
better than 1%. We then justify the choice of Zi — 100. 



D. Simulation box size 

As discussed by H08 and [8|, box size affects nonlin- 
ear power spectrum mainly through two effects. One 
is the sample variance. Smaller box size simulation suf- 
fers larger statistic fluctuations due to fewer independent 
modes. Our results are relatively insensitive to this sam- 
ple variance, because we take the ratio of the power spec- 
tra, which share more or less the same cosmic variance, 
thus the ratio will cancel much of this sample variance. 
The other is systematic errors induced by missing large- 
scale modes, which contribute to the tidal force. Ref- 
erence (see Fig. 9 of [Ij) shows that the systematic 
error is substantially small (< 1% at z = 0) even in box 
size L = llG/i-^Mpc out to fc IG/iMpc"^ We thus 
think that the 300/i~^Mpc box size is suitable. For the 
adopted box size, the largest available mode lies in the 
linear regime even at z = 0, allowing us to test the sim- 
ulation result against the linear evolution to check the 
modified GADGET-2 code. 

Based on the above tests, we justify that, with 512^ 
particles, a 300/i~^Mpc box size and initial redshift 
Zi = 100, N-body simulation based on GADGET-2 can 
help us reliably probe the matter power spectrum to 
k = 3/iMpc~^. The primary statistics that we inves- 
tigate in this paper, e, namely the ratio of power spectra 
of MG models and ACDM, can reach unity to an accu- 
racy of 1%. As we will find later, in MG models, even 



for those with moderate deviation from GR (e.g. 10% 
deviation), e can deviate from unity to 0(0. 1), an order 
of magnitude larger than the simulation error. We thus 
are confident that the resulting e is robust. 
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0.504 


1.000 


1.485 


1.972 


3.442 


Redshift 




0.313 


0.800 


1.269 


1.732 


3.142 






0.113 


0.600 


1.053 


1.500 


2.847 






0.00719 


0.500 


0.947 


1.383 


2.700 








0.400 


0.843 


1.271 


2.554 




N/A 




0.300 


0.741 


1.160 


2.415 








0.200 


0.639 


1.053 


2.281 






N/A 


0.150 


0.591 


1.002 


2.216 








0.100 


0.541 


0.951 


2.152 








0.050 


0.494 


0.901 


2.090 








0.000 


0.447 


0.850 


2.029 



TABLE I: The output redshifts for the simulation of each 
The baseline redshifts zs are that of C = 1 (ACDM). The 
corresponding redshifts Zi^ of ( ^ 1 are set up by D{zs,C, = 

1)=I>(2C,C). 



E. Checking Modified GADGET-2 

For the MG parameterization we adopt Eq. [2l we only 
need to do a minimal modification to the GADGET-2 
code [131 • The original code calculates the gravitational 
potential in GR. Multiplying it by a factor C, we obtain 
the potential %p in the MG models, which determines the 
acceleration of nonrelativistic particles. In GADGET-2, 
because the TreePM algorithm is adopted, the gravita- 
tional potential is explicitly split into a long-range part 
and a short-range part. We need to multiply both by the 
same factor 

We test the modified GADGET-2 code by comparing 
the simulated linear growth and the theoretically calcu- 
lated one. In the linear regime, the matter power spec- 
trum P(fc,z) oc D'^{z). The scale k = 0.025 /iMpc"^ is 
in the linear regime through all output redshifts, so we 
calculate the power spectrum at this scale and compare 
it to the theoretical prediction, given by Eq. [3] Fig. [3] 
shows the comparisons. The good agreement indicates 
that our modified GADGET-2 is correct. 
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T ' 1 ' 1 ' 1 — ] r 1 ' 1 ' 1 ' 1 — ] [ 1 ' 1 ' 1 ' r 




A'{k;{=1.0) a'(k;{=1.0) 

FIG. 4: The ratios of the nonhnear power spectra of MG models to the ones in ACDM when they have the same linear power. 
Results shown in different panels correspond to MG models of different ( values. Different colors in each panel show different 
output redshifts. The ones that reach larger A^(fc, C = 1) have lower redshifts. The two dotted black lines are in the 1% limit. 
In the bottom-right panel, we show the curves for e of different values of (shown in different line-styles). This comparison 
locates at the corresponding redshift with Zs = 1.2. [See the electronic edition of the Journal for a color version of this figure.] 



IV. SIMULATION RESULTS 

For the baseline — ACDM simulation, we choose 21 
snapshots, whose redshift zs is shown in Table H] We 
then run five MG simulations. We turn on the modi- 
fled gravity at z < zmg — 100. This is certainly not 
an unique choice, nor backed up by a solid argument. 
For example, we could turn on the modified gravity at a 
much later epoch. This will be the topic for future study. 
Naturally, ^ adopted in the simulations should cover the 
range allowed by present observations. Although there 
is no direct constraint for this particular type of MG in 
the literature, current constraints on MG (e.g. the pa- 
rameterization investigated by 15511 . the gravitational slip 
parameter w in [56j and r] in j57l | ) imply that, the cur- 
rent constraint on C reaches no better than ~ 10% accu- 
racy. This instructs us to adopt ( = 0.8, 0.9, 1.1, 1.2 for 
the simulations. We also run a simulation with C, — 1.5. 
Structure growth in the = 1.5 model is very likely too 
fast to fit existing observations (Fig. Nevertheless, 
we include this simulation since dramatic modifications 



in GR highlight signatures of MG in the nonlinear evo- 
lution, lead to better understanding of nonlinearity in 
MG models and help to improve the generality of fitting 
formulae based on these simulations [75j. 

The outputs of these simulations are chosen according 
to Eq. [21 such that the linear power spectrum of the 
given MG model at is identical to that of ACDM at 
Zs. The corresponding is shown in Table ID Since all 
simulations stop at z = and linear density growth in 
C < 1 universe is slower than that in ACDM (Fig. 
there will be no available z,^ for comparison in ^ < 1 
simulations. 

The main simulation results are shown in Fig. 2] As 
a reminder, the function e — 1 is the relative difference 
between the nonlinear matter power spectrum in the MG 
model and the corresponding one in the standard ACDM 
[Eq. [S]. By design, it scales out the difference in the 
linear density growth rate and thus highlights other fac- 
tors determining the nonlinear power spectrum. Further- 
more, due to this particular design, we are able to reduce 
the possible numerical artifacts and model the nonlin- 
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earity to 1% accuracy. To better show the effect of non- 
hnearity, we plot e against the nonhnear matter power 
spectrum variance A'^{k;zs,C = !)• To illustrate the 
evolution of e, we plot all e of the same C in the same 
plot. At the bottom-right panel of Fig. |4l we show a 
comparison between different values of C- The redshifts 
for different C of this comparison are selected when they 
have the same linear power spectrum as KCDM simula- 
tion at Zs = 1.20. 



A. Signatures of modified gravity in tlie nonlinear 
regime 

The results in Fig. |4]shows significant imprints (e 7^ 1) 
of modified gravity in the nonlinear evolution. Deviation 
of e from unity becomes stronger at smaller scales and 
lower redshifts where nonlinearity is stronger. We find 
that, by proper scaling of A^, curves of A^-e of fixed ( 
can fall upon each other. We will discuss this behavior 
further in the Appendix and show its application to de- 
velop a fitting formula of e, accurate to ~ 1%. Whether 
or not this behavior is generic will be investigated in fu- 
ture studies. 

There is another interesting behavior in the nonlinear 
evolution. The density growth in C > 1 cosmology is 
faster (Fig. However, after scaling out the linear 

growth, the normalized nonlinear evolution is actually 
slower (namely, e < 1 when C > Fig. |4]). When C < li 
the behavior is opposite (e > 1). 

The halo model may explain such behavior. The non- 
linear power spectrum can be decomposed into two terms 
fsi: 



A2 



NL 



(fc,z) = Ai(fc,z) 



°° dr) 

MS(k\M, z)b(M) dM 

dM 



Nr5^{k\M, z) 



dn 
dM" 



-dM 



(6) 



Here, S{k\M,z) is the Fourier transform of the density 
profile of a halo with mass M at redshift z, normalized so 
that S{k 0) — ?► 1. The halo abundance is dn/dM, and 
for convenience it is normalized such that / Mdn = 1. 
The first term on the right hand of the equation is the 
two-halo term, which dominates in the linear regime. The 
second term is the one-halo term, which dominates in the 
strongly nonlinear regime. Since S{k ^ 0) ^ 1 is inde- 
pendent of the halo density profile and J Mh{M)dn — 1, 
the two-halo term is not very sensitive to the halo den- 
sity profile. Since D{zq,Q = D{zs,( = 1), as required, 
the two-halo term in the MG model is (roughly) equal to 
the two-halo term in ACDM, where they are dominant. 
On the other hand, in the nonlinear regime, the one-halo 
term strongly depends on the halo density profile, which 
depends on the structure growth history. Since structure 
grows faster in C > 1 cosmology, we have z,^ > zg. Halos 
in this universe form in a background with higher mean 
density (oc (1 -I- z)~^). We then expect them to have a 



smaller concentration [69|. For the same mass, a smaller 
concentration means a smaller S{k\M,z) and a smaller 
contribution from the one-halo term. We then expect 
A^L(fc; z^X) < A^L(fc; zs, C = 1) and thus e(C > 1) < 1 
in the nonlinear regime. We defer this investigation to a 
forthcoming paper, where we will measure and compare 
the halo mass functions and profiles between MG models 
and ACDM. For the same reason, we expect e{( < 1) > 1 
in the nonlinear regime. Whether or not the halo model 
will lead to a satisfying description of e and thus the non- 
linear power spectrum in MG models is an interesting 
project for further investigation. 

In a word, large deviation of e from unity implies that 
there is valuable information of gravity encoded in the 
nonlinear regime, which is complementary to those en- 
coded in the linear matter power spectrum at the same 
epoch. It will be interesting to quantify how significantly 
this imprint of gravity in the nonlinear regime can im- 
prove cosmological tests of gravity. 



B. Implications on tlie applicability of some 
existing fitting formulae 

The particular definition of e allows us to address a key 
question in understanding the nonlinearity, is the non- 
linear power spectrum uniquely determined by the linear 
one at the same epoch? Equivalently, if the linear matter 
power spectra of two cosmologies are identical, will the 
corresponding nonlinear matter power spectra be identi- 
cal? 

The influential HKLM (Hamilton, Kumar, Lu and 
Matthews) procedure [58[ assumes so. It postulates the 
existence of an one-to-one mapping between the linear 
correlation function at a linear scale and the nonlin- 
ear correlation function at the corresponding nonlinear 
scale. Reference [11] found that this mapping depends 
on the slope of the linear power spectrum. Hence af- 
ter, the slope dependence has been explicitly incorpo- 
rated in several fltting formulae, including the popular 
PD96 fitting formula 



and the Smith et al. 2003 
halofit formula [47[. In these fitting formulae, the map- 
ping is expressed in Fourier space, of the functional form 
A^l(A;nl,2;) = ^(A|(fc,2;)). This mapping is nonlo- 
cal. For example, in PD96, A^l(^nl,-z) depends not 
only on A^(fcL,z) at the corresponding linear scale k^, 
but also on the effective power index ricff at some lin- 
ear scale, often chosen to be k^ or fcL/2. Further- 
more, the mapping has extra dependences on cosmol- 
ogy. In PD96, A^^^ik^i^,z) = H(Ai(fc, z), g(z)). The 
cosmological dependence g{z) has clear physical meaning, 
g{z) oc _D(z)(l -I- z) and is normalized to g{z — >■ 00) = 1. 
In the halofit, A^L(fcNL, z) = ^(A|(fc, z), flm), where flm 
explicitly enters several fitting parameters. For a com- 
prehensive review, refer to |47| . 

The HKLM procedure and its variations are successful 
in capturing the nonlinearities in CDM plus GR simula- 
tions. For this reason, they are often extended to predict 
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the nonlinear matter power spectrum in MG / dark energy 
models 153, [60l - l63| . The applicability of the resulting 
fitting formulae to MG models has been tested against 
several MG simulations [i^, [13, [fi^ ■ In general, there is 
reasonable agreement at the ~ 10% level. But discrepan- 
cies are also noticed (e.g. [13, [fi^). Our simulations, with 
improved simulation accuracy, improved power spectrum 
measurement method and specifically designed statistics, 
are able to identify the discrepancies at the 1% level. 

Our simulation results (Fig. 2]) show unambiguously 
that e 7^ 1 in the nonlinear regime. For a 10% devia- 
tion from GR (C = 1.1, or ( = 0.9), the resulting non- 
linear power spectra can differ by 20%-30% at i5 ~ 10 
(A^ ~ 100), to the corresponding one in the ACDM. The 
deviation becomes larger if the deviation of C from unity 
is larger. Quite obviously, the nonlinear matter power 
spectrum is not completely determined by the linear one 
at the same epoch. Reference [i^ demonstrated by the 
case of dynamical dark energy, that the structure growth 
history is also responsible for shaping the nonlinear mat- 
ter power spectrum. The MG models we investigate have 
a different structure growth history (e.g. different struc- 
ture growth rate), and this may explain the observed 
significant deviation of e from unity. 

Our simulation set up allows us to evaluate the appli- 
cability of using several existing fitting formulae to MG 
models even without directly testing them. Since all the 
simulations have identical linear power spectra and the 
present day matter density r2,„, the halofit would then 
predict e = 1. Our simulation result of e then implies 
that a ^ 20% error may occur if one uses the halofit to 
calculate A^l(C = 1-1) (^^nd A^l(C = 0.9)) at the over 
density 6 ~ 10. The application of PD96 to MG mod- 
els is a little bit tricky. In PD96, g{z) is the ratio of 
the linear density growth rate between the given CDM 
cosmology and the flm = 1 flat universe. However, this 
form of g{z), despite its clear physical meaning, does not 
apply to more general cases, such as the case of dynami- 
cal dark energy models [ll|. Furthermore, PD96 is based 
on the stable clustering hypothesis. N-body simulations 
show that this hypothesis is problematic [A7. . ,67] . Sim- 
ple extrapolation of PD96 to the MG models should be 
avoided too. 



V. DISCUSSION AND CONCLUSION 

In this paper, we modify the GADGET-2 TreePM code 
to run a set of simulations for parameterized modified 
gravity models. As the first paper in a series, we focus 
on the nonlinear power spectrum in MG models. We 
take several steps to improve/test the model and simu- 
lation accuracy. First, we adopt an advanced analysis 
method to improve the power spectrum measurement. 
We then test the impact of various mass and force res- 
olution, time step, initial redshift, and box size on the 
nonlinear power spectrum, and find suitable simulation 
specifications which meet our accuracy requirement. Fi- 



nally, we focus on a particular quantity e, the ratio be- 
tween the nonlinear power spectra between MG models 
and ACDM with the same power spectra at a large linear 
scale. This quantity can be measured to higher accuracy 
than the nonlinear power spectrum itself, since much of 
the sample variance and simulation artifacts are reduced 
in e. By construction, deviation of e from unity is a sig- 
nature of MG imprinted in the nonlinear evolution. It 
also means that the nonlinear matter power spectrum is 
not uniquely fixed by the linear one at the same redshift. 
It thus also represents the minimum systematical error 
induced by simply extrapolating some existing fitting for- 
mulae to these MG models. We find that, the deviation 
of e from unity can reach O(O.l) where the rms density 
fluctuation reaches 10, for MG models with a 10% devi- 
ation from GR. As an exercise toward a general fitting 
formula of this signature of MG, we develop a simple 
fitting formula of e, accurate to ~ 1%, working for the 
particular MG models that we investigate. 

Significant improvements are required to reach preci- 
sion modeling of the nonlinear matter power spectrum in 
more general MG models. In the next steps, we will run 
more simulations with larger box sizes, 1024^^ or more 
particles, various initial conditions, and various expan- 
sion histories. 

More importantly, we need to explore larger MG pa- 
rameter space. For example, we may need to vary zmg 
to see its influence. Furthermore, instead of taking as a 
step function with no scale dependence, we shall explore 
more complicated time dependent and scale dependent 

models. In a companion paper, we will explore the 
minimalist MG model (7- index, [3J|), which has been 
adopted by the Figure of Merit Science Working Group 
(FoMSWG) [zO] for forecasting. In this model, the linear 
density growth rate is given by /(a) = din D (a) /din a — 
n'^{a). Here, the growth index 7 is a constant, whose 
value is ~ 0.55 in ACDM. Stage-IV dark energy surveys 
can constrain this parameter with a rms error 0(0.01) 
(e.g. (Zli)- ^m{a) = rtoa'^ / {H^ / H^) is the matter den- 
sity at redshift z = 1/a — 1. One particular advantage 
of this parameterization is that, since Vlmio- — >■ 0) — >■ 1, 
even for a model with time-constant 7 can approach GR 
at high redshift. 

The C parameterization can incorporate this model. 
The corresponding <^ can be obtained from Eq. [31 

^_2 f + f{2 + H'a/H) + af' 

3 r2m(a) ' da 

It is interesting to see whether new features will arise 
in the nonlinear regime and how to extend the proposed 
fitting scheme for this MG model. 

Its is much harder to simulate realistic MG models 
such as f{R) and DGP, which have complicated envi- 
ronmental dependences. We hope that, studies on MG 
models without environmental dependence can provide 
useful templates to understand these MG models. 
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[72] defined here is equivalent to (in 3], up to a constant 
factor of AttG. 

[73] However, this does not mean that and the expansion 
rate uniquely fix the gravitational lensing statistics, since 
gravitational lensing is determined by the combination 
— -i/) cx Gofi {1 + l/rj)p = (^p{l + ri). Since this relation is 
algebraic, there is no extra simulations required to eval- 
uate the lensing statistics. The above arguments hold as 
long as Ge{!{k,z) and r]{k,z) are deterministic functions 
of scale k and redshift z. 

[74] This is not an arbitrary choice, as it appears to be. In gen- 
eral dark energy cosmology and MG cosmology, the ex- 
pansion history and the structure formation are indepen- 
dent. Since the structure formation is affected by quanti- 
ties such as the dark energy sound speed, the anisotropic 
stress, Gefi and rj (and () at sub-horizon scales do not 
affect the expansion. It is thus natural to fix the expan- 
sion rate matching the observations today and explore 
the possible difference in the structure growth. 

[75] We also run a simulation with ( = 0.5. We find that 
the structure growth in this model is too linear to be 
interesting. 



Appendix A: The fitting formula 

We demonstrate the feasibility to develop an accurate 
fitting formula for e, which quantifies the difference of 
the nonhnear evolution between the MG cosmology and 
ACDM. In combination with the existing ACDM fitting 
formulae, the nonlinear power spectrum in MG models 
can be predicted. We only use the results of C > 1 sim- 
ulations to develop the fitting formula. We reserve the 
C < 1 simulations to check the generality of our fitting 
formula. This fitting formula is by no mean applicable to 
general MG models. Nonetheless, we hope that it serves 
as an useful exercise toward more general fitting formula, 
which we will explore elsewhere. 



1. Developing the fitting formula 

As shown in Fig. 21 e is a function of the scale, redshift 
and C, 



cc = A2,L(fc,zs,C = l) . 



(Al) 



We will develop the fitting formula according to the fol- 
lowing steps. 

Step 1: By visually inspecting Fig. U] for each C, it 
looks feasible to move each lines horizontally such that 
they fall upon each other. This horizontal shift indeed 
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FIG. 5: Similar to Fig. IH but only shown for results with > 1. In addition, lines corresponding to different redshifts are 
shifted by a factor of s{z) as described in Eq. IA3I These lines are then fitted with three parameters A, B and C. The asterisks 
are the fitting points produced by the fitting recipe. 



works (Fig. [5]). Mathematically, this means that we can 
find a shift function s[zg, i^), such that 



u{x, zs, C) = u{y, C); y = xx s{zs, C) 



(A2) 



In the exercise, we fix ^ first. Then for each zg (or the 
corresponding z^), we find the suitable s, such that the 
corresponding curve, after the horizontal shift, overlaps 
with the one with zs = 0- To do so, we have required 
that s{zs = 0, C) = 1- We find that higher redshift curves 
should be shifted more leftward. This requires that (i) 
5(2:5) is positive, and (ii) s{zs) monotone decreases with 
respect to zs- These requirements help us to find the 
following fitting function for s: 



D{zs,C = l) 



i?(z5 = 0,C = l) 



(A3) 



where A{Q > 0. 

Step 2: The next step is to figure out a suitable form 
for u{y, C), namely, to fit those curves in Fig. [S] There are 
several guidelines. (1) u > 0, since both power spectra 
must be positive. (2) u{y,C = 1) = 1, by the definition. 
(3) e = u < 1 when C > 1 and e > 1 when C < 1 (Fig- 
2]). These behaviors motivate us to propose the following 
fitting function: 



=,(i-C)B(C)y^<^' 



(A4) 



As long as i?(C) > 0, all three conditions are satisfied. 

Step 3: For each (, we fit the simulation data and find 
the best fit A, B and C, whose values are shown in Fig. 
m We then need to find the suitable form to model the C 
dependence of these parameters. The following functions 
with the associated parameters provide a good fit: 

A{C} = e^""-'^'' , ao = 1.745 , 



B{0 = bo + hC\ 

c(C) = 0.573 . 



60,1=0.0429,0.133, (A5) 



Calculating the nonlinear power spectrum in 
MG models 



We summarize the procedure to calculate the nonlin- 
ear matter power spectrum A'^^(k, zq, C) using our fitting 
formula. 

• For the given redshift zq in the MG model, find 
the corresponding zs in ACDM through Eq. |4l 
The two corresponding power spectra at large lin- 
ear scale are then identical. 

• Calculate l\%^{k,zs,C, = 1). This can be done by 
using either direct ACDM simulations (as in our 
case) or existing fitting formulae such as the halofit. 



• Combining Eqs. [AH |A3l |A4] & |A5l we are then 
able to predict e(fc,z^,C). In combination with 
A^l(^j-^S',C = 1)1 we can then predict the non- 
linear matter power spectrum l\^^(k, zq^Q in the 
given MG model. 



3. Testing the fitting formula 

Becuase of the very limited simulations that we have, 
we are not able to perform comprehensive tests against 
the generality of the fitting formulae. However, we are 
indeed able to check it against our ^ < 1 simulations. 
Since we do not use these simulations to find the fitting 
parameters, these C < 1 simulations can provide an inde- 
pendent check against our fitting formula. 

In Fig. ini we show the performance of our fitting for- 
mula. The ratios between the simulated and predicted 
nonlinear power spectrum of MG models are shown using 
different line styles for different C, and redshifts as indi- 
cated. In general, our fitting formula is accurate to 1-2% 
in the range k < 3/i/Mpc. Although results for the cases 
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FIG. 6: Testing the accuracy of the fitting formula. We plot 
the power spectrum ratios as a function of k between the 
simulation results and our model predictions. Different ( with 
its checking redshift are shown in different colors and line 
styles. The vertical line is at fc = 3/iMpc^^, below which our 
simulation results are reliable, as shown in Sec. IIIII 



with C < 1 extrapolations of our fitting formula, their 
performances are as good as for the C > 1 cases. Such 
good performance demonstrates the applicability of our 
fitting formula. And we do not try to seek for possible 
slightly more accurate but much more complicated fitting 
formulae. 

We have shown that the proposed fitting formula pro- 
vides a good description of the nonlinear matter power 
spectrum, for the specific form of MG that we adopt. 
Is it applicable to other cases? We are not able to an- 
swer this question by the existing simulations. However, 
we still want to discuss a less general question, is it ap- 
plicable to the MG models with zmg 100? In the 
fitting formula, the only quantity dependent of zmg is 
D{zs,C = 1) = D{z(^,()- This dependence alone may 
not be sufficient. It is very likely that A, B, and C 
depend on zmg, too. This is a key issue for future in- 
vestigation. Nevertheless, we show that it is possible to 
develop a fitting formula for MG models by the above 
simple technique. It may also be applicable to more gen- 
eral MG models. This is again an interesting issue for 
further investigation. 



